mountie



An R package providing an efficient C++ implementation of the Reactive Multi-Particle Collisions (RMPC) algorithm.

Examples

library(mountie)

Diffusion-only system (MPC)

(system <- rdsystem_examples("diffusion"))
#> $network
#> NULL
#> 
#> $volume
#> [1] "S1"
#> # dims: 100 x 5 x 5
#> # h: 0.01
#> # states: 
#> # # A tibble: 2,500 x 4
#>       x     y     z    S1
#>   <int> <int> <int> <dbl>
#> 1     1     1     1     0
#> 2     2     1     1     2
#> 3     3     1     1     1
#> 4     4     1     1     0
#> 5     5     1     1     0
#> # โ€ฆ with 2,495 more rows
#> 
#> $D
#> S1 
#>  1 
#> 
#> $T
#> [1] 0.1
#> 
#> $tau
#> [1] 1e-04
#> 
#> $boundaries
#>       x            y          z         
#> lower "reflective" "periodic" "periodic"
#> upper "reflective" "periodic" "periodic"
#> 
#> $absorb_probs
#> NULL
#> 
#> $species_source
#> NULL
sol <- rmpc(system)
#> Starting RMPC simulation with the following parameters:
#>  - D:  
#>     S1 = 1
#>  - tau: 0.0001
#>  - T:   0.1
#> ....................................................................................................
sol$u[[1]]
#> # A tibble: 27,000 x 7
#>    Type         P.x     P.y     P.z    V.x    V.y    V.z
#>    <fct>      <dbl>   <dbl>   <dbl>  <dbl>  <dbl>  <dbl>
#>  1 S1-bath 0.00903  0.00259 0.00492  59.4   -10.4  -38.7
#>  2 S1-bath 0.0118   0.0115  0.00478 122.     33.3  -36.0
#>  3 S1-bath 0.0130   0.0477  0.0132   34.7   -33.1   76.1
#>  4 S1-bath 0.0291   0.0487  0.0349  248.    -17.4 -176. 
#>  5 S1-bath 0.00498  0.00673 0.0109   79.4    48.8   66.6
#>  6 S1-bath 0.00509  0.0347  0.0493   91.5  -190.   -31.2
#>  7 S1-bath 0.000214 0.0379  0.0478  -53.0  -163.   -91.3
#>  8 S1-bath 0.000769 0.0111  0.00577   3.59   43.1  -25.3
#>  9 S1-bath 0.0194   0.0124  0.0125  194.     53.9   33.0
#> 10 S1-bath 0.00364  0.0125  0.0381  106.     52.8 -168. 
#> # โ€ฆ with 26,990 more rows
px_plots <- lapply(1:length(sol$u), function(i) {
    position_plot(sol, system$volume, "S1", frame_index = i) +
        ggplot2::ylim(0, 2.5)
})
plot_indicies <- round(seq(1, length(px_plots), length.out = 3))
patchwork::wrap_plots(px_plots[plot_indicies], ncol = 1)

Gradient model with source

(system <- rdsystem_examples("morphogen"))
#> $network
#> # Reaction network: 1 reaction x 1 species
#>     Reactants    Products  Rate
#> 1           X -> 0            1
#> 
#> $volume
#> [1] "X"
#> # dims: 300 x 1 x 1
#> # h: 0.01
#> # states: 
#> # # A tibble: 300 x 4
#>       x     y     z     X
#>   <int> <int> <int> <dbl>
#> 1     1     1     1     0
#> 2     2     1     1     0
#> 3     3     1     1     0
#> 4     4     1     1     0
#> 5     5     1     1     0
#> # โ€ฆ with 295 more rows
#> 
#> $D
#>   X 
#> 0.5 
#> 
#> $T
#> [1] 3
#> 
#> $tau
#> [1] 1e-04
#> 
#> $boundaries
#>       x            y          z         
#> lower "reflective" "periodic" "periodic"
#> upper "reflective" "periodic" "periodic"
#> 
#> $absorb_probs
#> NULL
#> 
#> $species_source
#> $name
#> [1] "X"
#> 
#> $rate
#> [1] 6.25
#> 
#> $xrange
#> [1] 0 0
#> 
#> $yrange
#> [1] 0.00 0.01
#> 
#> $zrange
#> [1] 0.00 0.01
#> 
#> attr(,"class")
#> [1] "source"
sol <- rmpc(system)
#> Starting RMPC simulation with the following parameters:
#>  - D:  
#>     X = 0.5
#>  - tau: 0.0001
#>  - T:   3
#> ....................................................................................................
sol$u[[1]]
#> # A tibble: 6,000 x 7
#>    Type       P.x     P.y     P.z    V.x     V.y   V.z
#>    <fct>    <dbl>   <dbl>   <dbl>  <dbl>   <dbl> <dbl>
#>  1 X-bath 0.00568 0.00936 0.00809  -1.31   25.4  -90.8
#>  2 X-bath 0.0157  0.00974 0.00186  82.1   -69.6   11.5
#>  3 X-bath 0.00668 0.00321 0.00644   7.02   66.4  159. 
#>  4 X-bath 0.00860 0.00204 0.00480 104.     13.2  -36.6
#>  5 X-bath 0.00289 0.00198 0.00646  35.9  -101.    94.9
#>  6 X-bath 0.00527 0.00750 0.00122 -42.5    85.6  -42.6
#>  7 X-bath 0.00195 0.00794 0.00688  50.3    -4.68 -30.6
#>  8 X-bath 0.0201  0.00328 0.00148 139.   -118.   -46.9
#>  9 X-bath 0.0126  0.00655 0.00975  40.6   -25.9   44.0
#> 10 X-bath 0.00489 0.00180 0.00620 104.     58.6   21.6
#> # โ€ฆ with 5,990 more rows
px_plots <- lapply(1:length(sol$u), function(i) {
    position_plot(sol, system$volume, "X", frame_index = i) +
        ggplot2::ylim(0, 8)
})
plot_indicies <- round(seq(1, length(px_plots), length.out = 3))
patchwork::wrap_plots(px_plots[plot_indicies], ncol = 1)
#> Warning: Removed 300 rows containing missing values (geom_col).
#> Warning: Removed 245 rows containing missing values (geom_col).
#> Warning: Removed 232 rows containing missing values (geom_col).